
global c a_t b p e t n_t n_p x_t y_t ar nth0 pwr;

%*******************Set Geometry of the Domain***********************
% From command line 1) pdecirc(0,0,1); 2) pdecirc(0,0,0.7); 3) Subtract circles 
%4) clicka 'Partial Omega' 5) click 'Boundary' and Export geometry; 
% Export boundary condition to the matrix b 

%*********************You may select one of the following geometries***********

%pdecirc(0,0,1); %Disk
pdecirc(0,0,1), pdecirc(0,0,0.7); %annulus
%pdecirc(2,0,1), pdecirc(-1,0,0.5), pderect([-0.55 1.03 -0.2 0.2]); %Dumbbel
%pdetool('changemode',0);
%pdetool('removeb',[1,2,3,5,6,8,13,16]);
%pdecirc(2,0,1), pdecirc(-1,0,0.5), pdepoly([-0.55 0.25 1.03 1.03 -0.25 -0.55],[0.2 0.19 0.2 -0.2 -0.19 -0.2]); %Dumbbel
%pdeellip(2,0,.5,1,0),pdeellip(2,0,.1,.5,0),pdecirc(-1,0,0.5),pderect([-0.54,1.51,-0.2,0.2]);%Dumbbel
%pdepoly([-1 0 0 1 1 -1],[0 0 1 1 -1 -1]);                  %L-shaped polygon
%pdepoly([0 -0.5 0.5],[0.866025403 0 0]);     %unit triangle

%pderect([-.5 .5 -.5 .5]);

disp('Is the geometry correctly set? If no, please modify. If yes, press any key.');
pause;

[p,e,t,n_t,n_p,x_t,y_t,ar]=mesh_refine(g);

disp('Is the mesh O.K.? If no, please refine. If yes, press any key.');
pause;
